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ABSTRACT 

Aims. The aim of this paper is to study the efficiency of different approaches to interloper treatment in dynamical 
modelling of galaxy clusters. 

Methods. Using cosmologicaW-body simulation of standard ACDM model, we select 10 massive dark matter haloes 
and use their particles to emulate mock kinematic data in terms of projected galaxy positions and velocities as they 
would be measured by a distant observer. Taking advantage of the full 3D information available from the simulation, 
we select samples of interlopers defined with different criteria. The interlopers thus selected provide means to assess 
the efficiency of different interloper removal schemes found in the literature. 

Results. We study direct methods of interloper removal based on dynamical or statistical restrictions imposed on 
ranges of positions and velocities available to cluster members. In determining these ranges, we use either the velocity 
dispersion criterion or a maximum velocity profile. We also generalize the common approaches taking into account 
both the position and velocity information. Another criterion is based on the dependence of the commonly used virial 
mass and projected mass estimators on the presence of interlopers. We find that the direct methods exclude on average 
60-70 percent of unbound particles producing a sample with contamination as low as 2-4 percent. Next, we consider 
indirect methods of interloper treatment which are applied to the data stacked from many objects. In these approaches, 
interlopers are treated in a statistical way as a uniform background which modifies the distribution of cluster members. 
Using a Bayesian approach, we reproduce the properties of composite clusters and estimate the probability of finding 
an interloper as a function of distance from the object centre. 

Key words, galaxies: clusters: general - galaxies: kinematics and dynamics - cosmology: dark matter 
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■ 1. Introduction with velocities larger than Suios. This simple approach is 

k> , still widely used today. With enough galaxies in a sam- 

. The modelling of galaxy kinematics in clusters remains one pie, one can take into account the dependence of ctios on 

d ' of the major tools in determining their properties, in partic- the projected distance from the cluster centre R and per- 

■ ■ ■ ' ular their mass distribution and dark matter content. Due form the rejection procedure in bins with different trios or 

to projection effects, any cluster kinematic data sample in- fit a simple solution of the Jeans equation to the measured 

evitably contains galaxies that are not bound to the cluster line-of-sight velocity dispersion profile, aios{R), and reject 

and therefore are not good tracers of its gravitational po- galaxies outside the 3crios(-R) curves (Lokas et al. 2006). 

tential. We will call these galaxies interlopers. An essential pg^g^ et al. (1990) discussed another method relying on 

step m dynamical modelhng of clusters by any method is iterative removal of galaxies whose absence in the sample 

therefore to remove such interlopers from the samples or ^^^^^g ^jgg^g^ change in the mass estimator. Zabludoff 

take their presence into account statistically Velocity mfor- (^ggg)^ Katgert et al. (1996) and Fadda et al. (1996) 

mation can be use_d to remove obvious interlopers that are advertised the use of gaps in the velocity distribution as 

thousands of km s off the mean cluster velocity, but there ^ separate interlopers from real cluster members, 

remam numerous interlopers that he m a similar general oiaferio & Geller (1997) and Diaferio (1999) proposed the 

velocity range as the cluster members. Some hints can be caustics where the projected distribution function is 

provided by studying the photometric properties of galax- sufficiently low to separate cluster members from the sur- 

les or restricting the samples to elliptical galaxies but these mounding medium. Prada et al. (2003) discussed the solu- 

approaches usually do not solve the problem completely. ^ion to the problem based on the use of escape velocities. 

It has long been recognized that the line-of-sight veloc- The first methods that combine the information on the po- 

ity distribution of galaxies in clusters is close to a Gaussian, sition and velocity of a galaxy were proposed by den Hartog 

The first attempts to design a scheme to remove the inter- & Katgert (1996) and Fadda et al. (1996). All these meth- 

lopers were based on this property. Yahil & Vidal (1977) ods aim at cleaning the galaxy sample from non-members 

proposed to calculate the line-of-sight velocity dispersion before attempting the proper dynamical analysis of the 

of the galaxy sample, trios, and iteratively remove outliers cluster; we call them direct methods of interloper removal. 
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A completely new approach to interloper treatment was 
pioneered by van der Marel et al. (2000) where, for the first 
time, the interlopers were not identified and removed from 
the sample, but their presence was taken into account sta- 
tistically by appropriate modification of the distribution 
function of the galaxies. A similar approach was also con- 
sidered by Mahdavi & Geller (2004) with more realistic as- 
sumptions concerning the distribution of interlopers. Prada 
et al. (2003) studied the distribution of satellites around 
giant galaxies by fitting to the projected velocity distri- 
bution the sum of a Gaussian and a uniform distribution 
taking care of the background. We will refer to this type of 
methods as indirect. It should be noted that these meth- 
ods are mainly applicable to composite clusters, i.e. data 
sets created by combining kinematic data from many ob- 
jects because only then the samples are numerous enough 
to provide useful constraints on the interloper fraction. 

The different methods of interloper treatment found in 
the literature are difHcult to compare. Each one of them has 
a different set of underlying assumptions. They also differ 
by the amount of parameters that have to be put in by 
hand. Most of the methods are iterative and some may not 
converge. The ultimate comparison between the methods 
can only be performed by resorting to A^-body simulations 
where full 3D information is available and true interlopers 
can be identified. Such tests have been already attempted 
(e.g. by Pcrca et al. 1990; den Hartog & Katgcrt 1996; 
Diaferio 1999). In particular, van Haarlem et al. (1997) 
compared the methods of den Hartog & Katgert (1996) 
and Yahil & Vidal (1977) in terms of the quality of re- 
production of the real velocity dispersion. However, more 
systematic study of different procedures is still needed and 
this is the aim of the present paper. We implement and gen- 
eralize different prescriptions for interloper removal found 
in the literature and apply them to mock kinematic data 
created from the simulation. Our goal is to measure the ef- 
ficiency of the different methods by measuring fractions of 
interlopers they remove. 

Our choice of methods will of course be arbitrary. We 
tried to focus on those easiest to implement, most widely 
used in the literature and with the smallest number of pre- 
selected parameters so that they are applicable not only to 
galaxy clusters but to all astronomical systems where kine- 
matic measurements of discrete tracer can be made (e.g. 
dwarf spheroidal galaxies). In the near future we plan to 
apply the methods discussed here to nearby clusters from 
the WINGS survey (Fasano et al. 2006) where about 300 
redshifts per cluster will be available. 

The problem of the treatment of interlopers is directly 
related to the problem of the mass estimation in gravita- 
tionally bound objects. We will demonstrate in section 2 
that using contaminated kinematic samples can lead to se- 
rious errors in the estimated mass. In addition, several of 
the interloper removal schemes we discuss make use of some 
crude mass estimators. However, the purpose of this work 
is not to provide the best method for mass estimation in 
galaxy clusters. Instead, we focus on a much narrower is- 
sue of how to obtain a clean sample of cluster galaxies free 
of interlopers before attempting a further analysis of the 
mass distribution in the cluster. This final analysis can be 
performed via a number of methods e.g. fitting velocity 
dispersion profile assuming isotropic orbits (e.g. Biviano & 
Girardi 2003), fitting velocity dispersion and kurtosis for 
arbitrary constant anisotropy (Lokas et al. 2006) etc. The 



final outcome of these procedures will depend on their spc;- 
cific properties and on the properties of objects to which 
they are applied (e.g. whether they are spherically sym- 
metric, depart from equilibrium, how well they are sam- 
pled etc.). For example, Sanchis et al. (2004) and Lokas et 
al. (2006) apphed the dispersion-l-kurtosis fitting method 
to simulated clusters (after removal of interlopers) and dis- 
cussed how well the main properties of the clusters (includ- 
ing the mass) are reproduced. 

For the purpose of this study we used a present-day 
output of a pure dark matter, medium-resolution cosmo- 
logical Af-body simulation in which cluster-size haloes can 
be identified. Taking advantage of the fact that the distri- 
bution of galaxies in clusters is similar to mass distribution 
in simulated dark matter haloes, i.e. both are cuspy and 
can be approximated by the NFW (Navarro et al. 1997) 
profile (e.g. Carlberg et al. 1997; Lokas & Mamon 2003; 
Biviano & Girardi 2003) we assumed that the galaxies can 
be approximated by just dark matter particles. Although it 
would be worthwhile to test the methods on a set of higher 
resolution simulations where galaxies or subhaloes can be 
identified, the distributions of subhaloes both in space and 
velocity are known to be biased with respect to those of 
dark matter particles (Diemand et al. 2004). On the other 
hand, Faltenbacher & Diemand (2006) have recently shown 
that subhaloes with sufficiently high mass corresponding 
to galaxies have distributions much less biased and very 
similar to those of dark matter particles, which makes the 
effort of using subhaloes questionable. Nevertheless, a pos- 
sibility to assign stellar populations to subhaloes can con- 
siderably improve their usefulness in the analysis. As shown 
by Biviano et al. (2006), subhaloes with old stellar popu- 
lations are more concentrated around their mother haloes 
so by selecting them one can reduce the contamination by 
interlopers. 

The paper is organized as follows. In section 2 we de- 
scribe the way to create mock data sets from the simu- 
lations, introduce different types of interlopers and discuss 
how the presence of interlopers can affect the inferred prop- 
erties of a galaxy cluster. Section 3 is devoted to direct 
methods of interloper removal. We first discuss the dynam- 
ical approach where the maximum velocity available to a 
member galaxy is estimated using some assumptions about 
the cluster mass profile. Next we study the statistical ap- 
proach in its most commonly used forms which we then 
generalize by considering the distribution of galaxies in pro- 
jected phase space. We also discuss the efficiency of different 
mass estimators in identifying interlopers. Section 4 is de- 
voted to indirect methods of interloper treatment and the 
discussion follows in section 5. 

2. Interlopers on velocity diagrams of simulated 
haloes 

In this work we used an A^-body cosmological simulation of 
standard ACDM model described in Wojtak et al. (2005). 
The simulation was performed using a version of the ART 
(Adaptive Refinement Tree) code (Kravtsov et al. 1997) 
in a box of size 150 Mpc with parameters h — 0.7, 
r^M = 0.3, JIa = 0.7 and as = 0.9. From the whole sample 
of dark matter haloes formed in the final simulation output 
{z = 0) we choose 10 massive (lO-^^-lO^^M©) and possibly 
relaxed ones i.e. without any obvious signatures of ongoing 
major mergers. They are listed and described in Wojtak et 
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Fig. 1. Velocity diagrams of halo 6 out to i? = in the mass centre rest frame of reference seen in different projections 
and with different types of interlopers. FiUed and empty circles indicate halo particles and interlopers respectively. In 
the top left corner of each panel we mark the projection axis and the criterion for interloper identification (> r^, > 2rv 
or > Vcsc for particles beyond r^, 2r„ and unbound particles respectively). In the bottom left corner we give numbers of 
halo particles and interlopers which are seen on a given velocity diagram. 



al. (2005). All of them are characterized by mildly radial 
particle orbits and their density profiles are well fitted up 
to the virial radius r„ by the NFW formula 

P _ Acc^ff(c) 
Pcfi 3(r/r,)(l + c(rA-„))2' ^' 

where 5 (c) = [ln(l + c) — c/(l + c)]^^, c is the concentration 
parameter, is the critical density at present and Ac 
is a parameter defining virial mass in terms of overdensity 
with respect to the critical density. We assume Ac = 101.9 
which is the value valid for the concordance ACDM model 
with rim = 0.3 and I^a = 0.7 (Lokas & Hoffman 2001). The 
mean value of the concentration parameter averaged over 
all 10 haloes is equal to 7.2. 

In order to emulate kinematic data for a galaxy cluster 
embedded in a given dark mater halo we place an imagi- 
nary observer at the distance of D=100 Mpc from the halo 
centre (going from the simulation comoving coordinates to 
the observer's redshift space) so that the receding veloc- 
ity of a halo mass centre observed by him is around 7000 
km s~^. Approximating the conical shape of the observa- 
tion beam with a cylinder {D » r„), we project position 
vectors of simulation particles onto the plane perpendicular 
to the line of sight and their velocities with respect to the 
observer onto his line of sight. Assuming that some of the 



simulation particles represent galaxies, we randomly select 
300 particles from the inside of the observation cylinder 
with projected radius R = r^, where the virial radius r„ 
is found in 3D analysis. Additionally, we restrict our selec- 
tion to particles with velocities from the range ± 4000 km 
s~^ with respect to the velocity of a halo mass centre. This 
choice of velocity cut-off, corresponding to at least 4(Tios for 
cluster-size objects, guarantees that we do not exclude any 
cluster galaxies with high peculiar velocities. 

We place the cylinder of observation along the main axes 
of the simulation box so the orientation of the haloes (which 
have triaxial shapes) with respect to the line of sight should 
be random. Finally for each of the 10 haloes we obtain 
three sets of projected galaxy positions and velocities from 
observations along x, y and z axis of the simulation box. 
We will refer to these sets of data as velocity diagrams. Each 
velocity diagram includes both particles from the inside of a 
given halo (we call them simply halo particles) and particles 
from the outside of a halo which are seen because of the 
projection effects (we call them interlopers). 

First we identify the true interlopers in our data using 
the full 3D information about positions and velocities of 
simulation particles. Obviously one can think about inter- 
lopers as particles which are beyond r„ since they are from 
the outside of virialized region and they are not used in the 
estimation of density profile. We find that on average 24 
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percent of particles on our velocity diagrams have r > r„. 
This criterion of interloper identification, however, seems 
to be too restrictive in many cases since the object may 
possess a virialized region or at least a well defined density 
profile extending up to radii beyond defined by density 
contrast parameter Ac (e.g. Klypin et al. 2003; Wojtak et 
al. 2005; Prada et al. 2006). Besides, halo shapes are not 
spherical and imposing this kind of symmetry we can lose 
particles that are inside the real virialized region of the 
halo. We find that almost half of particles beyond Vy on 
our velocity diagrams reside below ~ 2r„ and are bound to 
their halo (the fraction of unbound particles is negligible 
at 2r^). These considerations suggest that they could also 
be treated as good tracers of the halo potential. We have 
therefore decided to consider two more conservative criteria 
of interloper selection: particles beyond 2r„ and unbound 
particles (with velocity greater than the escape velocity). 
Average contribution of these groups to the particles on the 
velocity diagrams is 13 and 8 percent respectively. Fig. [T] 
shows a set of velocity diagrams for halo 6 in different pro- 
jections (rows) and with different criteria of interloper se- 
lection (columns) , where filled and empty circles correspond 
to halo particles and interlopers respectively. 

To illustrate how interlopers affect the results of dy- 
namical analysis we fit for simplicity an isotropic solution 
of the Jeans equation to a velocity dispersion profile mea- 
sured for one of the simulated velocity diagrams (see Lokas 
& Mamon 2001, 2003 and Lokas et al. 2006 for details of the 
Jeans formalism and the fitting procedure) . The dispersion 
profile is measured in radial bins for a whole sample (300 
particles) and for three subsamples of particles cleaned of 
three types of interlopers introduced above. Fig. [5] shows a 
typical velocity diagram generated from our halo 4 (with 
filled and empty circles as bound and unbound particles re- 
spectively), dispersion profiles for the four mentioned sub- 
samples of particles and the corresponding results of the 
fitting procedure aimed at estimating the virial mass 
and concentration c of the halo, where parameter values 
found in 3D analysis are marked with a cross. All lines cor- 
responding to the same particle sample are drawn with the 
lines of the same type. 

Although the results shown in Fig. [2] concern just a sin- 
gle case of a velocity diagram, they illustrate well the gen- 
eral feature of bias caused by interlopers. First, note that 
the velocity dispersion is overestimated mainly in the outer 
part of the velocity diagram and this is caused mostly by 
unbound particles (since all dispersion profiles calculated 
for the data cleaned of interlopers of three different types, 
which include at least all unbound particles, are almost the 
same). Second, all three corrected dispersion profiles infer 
fitting results which are very similar to each other and in- 
clude the true parameter values inside la confidence level 
contour. Adding unbound particles to the analysis shifts My 
towards higher masses (which is due to the overestimated 
velocity dispersion) and forces the concentration parame- 
ter to lower values (which is due to the rising dispersion 
profile) . 




1 1.5 

R [Mpc] 




Fig. 2. The top panel shows the velocity diagram of halo 
4 in projection along z axis. Filled and empty circles indi- 
cate particles bound and unbound to this halo respectively. 
Middle and bottom panels show respectively the disper- 
sion profiles and the results of the fitting procedure in the 
form of 68.3, 95.4 and 99.73 percent probability contours in 
the My — c parameter plane assuming isotropic orbits. The 
different types of lines correspond to different subsamples 
of particles used to calculate the dispersion profile: dotted 
lines are for the whole sample with interlopers included, 
solid ones for bound particles, dashed ones for particles be- 
low 2r u and dotted-dashcd lines for particles below r„ . The 
cross marks the concentration parameter and virial mass 
found in 3D analysis of the mass distribution in the halo. 



3. Direct methods of interloper removal 

3.1. Overview 

In this section we study methods which allow us to remove 
a significant fraction of interlopers using some criteria. First 



we consider restrictions on the positions of halo particles on 
the velocity diagram. Given the maximum velocity available 
for halo particles (dynamical approach) or a distribution of 
halo particles on the velocity diagram (statistical approach) 
we impose boundaries on the area of the velocity diagram 
likely occupied by halo particles. Interlopers are then iden- 
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tified as particles from the outside of this area. Then we 
consider a criterion based on the way interlopers affect dif- 
ferent mass estimators. In all these approaches the proce- 
dure of interloper removal is iterative. In each step new 
boundaries of the area occupied by halo particles or mass 
estimators are determined from the data partially cleaned 
of interlopers in the previous steps and the next group of 
interlopers is removed. All methods are supposed to con- 
verge after a few iterations when no more interlopers are 
identified. 

Knowing which particles on the velocity diagrams are 
real interlopers (belonging to any of the samples defined 
in the previous section) we are able to study the efficiency 
of different methods aimed at eliminating interlopers from 
velocity diagrams by comparing lists of interlopers found 
by these methods with those identified in 3D analysis. To 
quantify these results we introduce three parameters: a frac- 
tion of identified interlopers ft , a fraction of halo particles 
(galaxies) which were taken for interlopers by mistake fg 
and a fraction of non-identified interlopers remaining in the 
final sample of halo members fh- For an ideal method of 
interloper removal we would have all interlopers identified 
correctly, i.e. fi = 1, fg = and fh — 0. In order to judge 
the performance of different schemes of interloper removal 
we calculate the mean values of the parameters and their 
dispersions averaging over the whole set of velocity dia- 
grams. It should be kept in mind, however, that the values 
of these parameters will depend on the initial velocity cut- 
off used to select the data (allowing wider velocity range 
we would obtain higher values of fi). The important point 
is that the relative efficiency of different methods of inter- 
loper removal should not depend on this velocity cut-off. 
We address this issue further in the last section. 

3.2. Dynamical approach 

In this approach, we identify an interloper as a particle at 
a given projected radius R whose velocity exceeds a maxi- 
mum velocity available for halo particles at this radius. The 
main problem of this method lies in the choice of proper 
maximum velocity profiles. Let us consider two characteris- 
tic velocities: the circular velocity Vdi and the infall velocity 
Vini given respectively by 



Vcir = VGM{r)/r 



(2) 
(3) 

Another quantity of interest would be the escape veloc- 
ity Vcsc = ^2|(f>(r)|, however, as discussed by Prada et 
al. (2003) it does not lead to any useful criterion for in- 
terloper removal. The interpretation of the infall velocity 
is as follows. Assuming circular orbits of a given set of 
particles one can obviously recover the relation between 
potential and kinetic energy {U and T respectively) pos- 
tulated by the virial theorem: 2Tcir = — t/cir- The infall 
velocity Winf is simply an upper limit to the particles' ve- 
locities for which the virial theorem equation is violated, 
Tinf = -U (den Hartog & Katgert 1996; Beers et al. 1982). 
This limit originates from the requirement that a given par- 
ticle is bound to the halo {U + T < 0). Note that this 
condition provides a stronger restriction on the maximum 
velocity than the general formula for the escape velocity 
since Wcir < •\/|<I'(r)|. The equality Ucsc — Vinf would only 
occur if the density distribution dropped to zero at r since 




Fig. 3. Maximum velocity profiles. The dark gray area in 
the background indicates maximum velocity reached by 
bound particles. Medium and light gray strips correspond 
respectively to the formula ([4]) and ([5]) with the mass pro- 
file given by the mass estimator ([5]). Widths of the shaded 
areas are given by the dispersions following from averaging 
over the whole sample of velocity diagrams. 



then we would have |$(r)| = GM{r)/r. The velocity Winf 
can therefore be viewed as an escape velocity from the mass 
interior to r. 

Following den Hartog & Katgert (1996) we introduce 
two formulae for the maximum velocity profile. First, as- 
suming that the direction of particle velocity in the limit 
determined by Winf has any orientation, the maximum ve- 
locity at a given projected radius R is given by 



maxi^{ui„f}, 



(4) 



where max^j is a maximum along the line of sight at the 
distance R from the halo centre. A second, more restrictive 
criterion which gives more accurate limits at high i? ~ r„ 
can be obtained from 



maxi^juinf cos6', Wcir sin 6*}, 



(5) 



where 9 is the angle between position vector of the particle 
with respect to the halo centre and the line of sight. With 
this formula, we assume a special kinematic model which 
allows particles to fall onto the halo centre with velocity Winf 
or to move in a tangential direction with circular velocity 

Weir- 

To complete the above prescription for the maximum 
velocity profiles one needs to specify the mass profile. As 
proposed by den Hartog & Katgert (1996), we use the mass 
estimator Mvt derived from the virial theorem (Limber 
& Mathews 1960; BahcaU & Tremaine 1981; Heisler et al. 
1985) 



MvT{r = i?n 



37riV Y^iivi-vf 
2G I]j<jl/i?ij- 



(6) 



where iV is a number of galaxies enclosed on the sky by 
a circle with radius i?max, Vi is the velocity of the i-th 
galaxy and Rij is a projected distance between i-th and j- 
th galaxy. This formula is valid for spherical systems with 
arbitrary anisotropy. The mass profile can be simply ob- 
tained as M{r) « MvT{Ri < r < Ri+i), where Ri is the 
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Fig. 4. Illustration of successive steps in Wmax(l) method 
of interloper removal for halo 6 in projection along the y 
axis. The method uses the maximum velocity profile ^ 
and mass profile ([5]) . The top and bottom panels show suc- 
cessive mass profiles and maximum velocity profiles sep- 
arating interlopers from halo particles respectively. Filled 
and empty circles mark particles bound and unbound to 
the halo. Numbers indicate successive steps of the proce- 
dure which are described by the lines of different types. 
Final maximum velocity profiles (number 4) are drawn with 
dashed-dotted lines. 



sequence of projected radii of particles (galaxies) in the in- 
creasing order. The virial theorem applies to a whole system 
and otherwise one needs to add a surface term (e.g. The & 
White 1986). Recently Biviano et al. (2006) estimated its 
value for simulated clusters within an aperture of 1.5 
Mpc. However, since our purpose is not to estimate accu- 
rately the mass but to design a procedure for interloper 
removal, the formula ([6]) is sufficient. 

Using subsamples of bound particles gathered on all 
mock velocity diagrams we calculate the maximum veloc- 
ity profiles given by ([4]) and (O with the mass profile de- 
termined by (III). The results, expressed in units of circu- 
lar velocity at the virial radius Vy, are shown in Fig. [3] 
as medium and light gray strips respectively, whereas the 
dark gray profile seen in the background of this plot in- 
dicates the average maximum velocity reached by bound 
particles (including bound particles beyond r^) on any ve- 
locity diagram. Widths of all three areas are given by the 
dispersions resulting from averaging the maximum velocity 
profiles over the whole sample of velocity diagrams. It is 
clear that formula ^ is expected to work best in remov- 



Table 1. Results of different procedures of interloper re- 
moval in terms of the fraction of removed interlopers fi, 
the fraction of halo particles incorrectly identified as in- 
terlopers fg and the fraction of non-identified interlopers 
remaining in the final sample of halo members fh- The ta- 
ble lists both the mean values (/) and the dispersions cr/ of 
the parameters obtained in the analysis of 30 mock veloc- 
ity diagrams. For all seven methods of interloper removal 
considered in section 3 results are quantified for three dif- 
ferent definitions of interlopers, particles beyond r„ or 2ry 
and unbound particles, as marked in the second column 
(interloper type - i/t) by r^, 2r^ and Vcsc respectively. 

method i/t (/,) a/, {fg) cr/^ {fn) Of^ 

"JwIT] ~v 23 17 n O 1977 8?^ 

2r„ 48 26 1.1 1.4 7.7 7.8 

W 73 23 1.0 1.4 2.4 4.0 

t;max(2) ~„ 13 n 00 00 21:9 

2r„ 30 22 0.0 0.0 10.2 9.0 

tJcsc 48 29 0.0 0.0 5.2 6.5 

3crios(5) 17 9 09 O 20 M~ 

2r„ 37 21 0.8 1.2 9.4 8.3 

t;csc 58 27 0.7 1.0 4.3 5.3 

3aios(10) rl 19 10 O 2A 2T:0 8J~ 

2r„ 40 22 1.5 2.4 9.1 8.3 

W 63 27 1.4 2.2 3.9 5.1 

3aios(ii) ~v 19 17 02 05 204 

2r„ 40 25 0.2 0.4 8.6 8.4 

TJcBc 61 28 0.2 0.3 3.4 4.4 

vihn r„ 19 17 0.3 0.5 20.4 9.2 

2r„ 41 25 0.3 0.5 8.5 8.3 

^cBc 62 27 0.3 0.4 3.4 4.4 

Mp/MvT 18 7 1.2 1.2 21.1 8.9 

2r„ 40 22 1.2 1.3 9.3 8.3 

TJcsc 65 26 1.2 1.2 4.1 5.8 



ing interlopers since its profile (light gray) coincides almost 
exactly with the maximum velocity reached by bound par- 
ticles (dark gray). On the other hand, the profile generated 
by formula (jl]) seems too conservative to be useful. 

Fig. m illustrates successive steps of interloper removal 
with the maximum velocity ^ (hereafter t)max(l) method) 
for one of mock velocity diagrams with rather large num- 
ber of unbound particles. The top and bottom panels show 
mass profiles and maximum velocity profiles separating in- 
terlopers from halo particles on the velocity diagram for 
successive iterations of this method marked with numbers. 
The final virial mass given simply by the value of estimator 
MvT for i?max = fy is cqual to 8.35 x IQ^'^Mq which is a few 
times lower than for the total contaminated sample (first 
iteration in Fig. 2]) and reasonably close to the real value of 
the virial mass found in 3D analysis, 5.35 x IO^^^Mq. Note 
that the mass estimator Mvt is known to overestimate 
the true mass due to the neglect of the surface term (see 
Biviano et al. 2006) and a more reliable final estimate can 
in general be obtained by fitting velocity moments (Sanchis 
et al. 2004; Lokas et al. 2006). 

With this method, on average 73 percent of unbound 
particles are identified and removed from a sample and only 
around 1 percent of bound particles are taken for interlopers 
and lost from the velocity diagram so that the final samples 
include only around 2-3 percent of unbound particles (see 
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Table[T]for details). Note that the fraction of removed inter- 
lopers fi is limited in principle to values lower than about 
75 percent because roughly 1 /4 of unbound particles within 
the observation cylinder with velocity cut-off 4000 km s~^ 
are within the envelope of bound velocities and therefore 
inaccessible for direct methods of interloper identification. 
Since fi = 73 percent available in this approach is very close 
to the expected maximum, the method presented above is 
possibly the most effective. As expected from Fig. [3l the 
method of interloper removal based on profile ^ is too 
conservative and on average identifies much less interlopers 
than the previous one (see Wmax(2) method in Table [J). 

3.3. Statistical approach 

The idea of this approach is to use the information about 
the distribution of halo particles on the velocity diagram to 
distinguish between the probable halo particles and inter- 
lopers. The first scheme along these lines was introduced 
by Yahil & Vidal (1977) who proposed to identify inter- 
lopers as galaxies with velocities from the outside of the 
range ±3crios around the mean cluster velocity, where ctios 
is the projected velocity dispersion of galaxies in the cluster 
given by the standard unbiased estimator. In this formula- 
tion the method is model-independent so that the data are 
self-verified as far as the interloper removal is concerned. 

It is easy to generalize the above prescription to the 
case of data gathered in n radial bins so that Sfiios proce- 
dure could be applied in each bin independently in the way 
proposed by Yahil & Vidal. This modification allows us to 
take into account dependence of the velocity dispersion on 
R. However, increasing the number of bins we let the dis- 
persion in the outer part of velocity diagram be much more 
overestimated by interlopers. A way to overcome this prob- 
lem is to use subsequently different numbers of bins. The 
dispersion in wide bins (when there is a small number of 
bins) is less biased by interlopers so in this case we remove 
interlopers efficiently. On the other hand, using narrow bins 
(when there is a larger number of bins) we measure the dis- 
persion locally taking into account the dependence of ctios 
on R. 

In each step of this method we use the following esti- 
mators of mean velocity and velocity dispersion 



m - 2 



(7) 
(8) 



where j is a number of bin, m is a number of datapoints 
per bin and Vij is the sequence of velocities in the i-th bin 
with the most outlying from the mean value in the last 
position so that following the prescription of Yahil & Vidal 
we do not take into account these velocities in estimating 
the dispersion. For each number of bins changing in a given 
step from ri,nin to ri,nax we remove particles with velocities 
from the outside of the range ±3(Tios,i around Vi : | Wi j — 5i | > 
ScTios.i- The procedure converges after a few steps when no 
more interlopers are found in any bin. 

We find that in order to remove even strongly clustered 
groups of interlopers (like the ones on the velocity diagram 
of halo 6 in projection along y axis seen in the middle left 
panel of Fig. [6]) it is necessary to fix nmin = 1 which corre- 
sponds to the original approach of Yahil & Vidal (1977). To 




Fig. 5. Dependence of the results of the interloper removal 
method based on the analysis of the velocity dispersion of 
binned data on the maximum number of bins rtmax- The left 
panel shows the fraction of identified unbound particles fi 
and the right one the fraction of bound particles taken for 
interlopers by mistake fg. Widths of both profiles (dark 
gray) are given by the dispersions following from averaging 
over the sample of velocity diagrams. Light gray strip in 
the background indicates the best range of Jimax- 



fix a maximum number of bins nmax we consider the depen- 
dence of fi and fg on different choices of the value of rimax- 
The results are shown in Fig.[5]in the form of dark gray pro- 
files. Using them we can easily find the values of rimax for 
which the procedure gives possibly high fi and low fg. This 
range of rimax is marked with a light gray strip in the back- 
ground of the plot. Applying its lower limit (rimax — 5), for 
which fi profile begins inclining towards rimax axis, leads 
to slightly conservative method of interloper removal with 
average rate of unbound particle identification fi — 58 per- 
cent (see 3crios(5) method in Table [1]). In the upper limit of 
this range (nmax — 10), when fg starts increasing rapidly, 
algorithm is a bit more restrictive and allows to remove 
on average 63 percent of unbound particles with the rate 
of misidentification fg comparable to the result of i'max(l) 
method (see 3crios(10) method in Tabled] for more details). 

Recently Lokas et al. (2006) generalized the ±3(Tios rule 
of interloper identification to a continuous velocity dis- 
persion profile: ±3ctios(^), where aiosiR) is the projected 
isotropic solution to Jeans equation parametrized by 
and c (Binney & Mamon 1982; see also Prugniel & Simien 
1997; Lokas & Mamon 2001; Mamon & Lokas 2005) fitted 
to the binned data. Assumption of isotropic orbits allows us 
to break the degeneracy between c and the anisotropy and 
to trace accurately the shape of the velocity dispersion pro- 
file with the c parameter only. We find that in some cases 
of velocity diagrams with strong interloper contamination 
this procedure stops too early because of the overestima- 
tion of the velocity dispersion proffie. To fix this problem 
we propose to fit aiosiR) to an incomplete dispersion pro- 
file after rejecting a few outer data points which are most 
contaminated. In our case we proceeded with the fitting 
for at least 6 data points to the maximum of 10 (always 
with 30 particles per bin) and then used the mean values 
of Aft, and c obtained for fc = 6 — 10 data points. The 
mean values were weighted with the goodness of fit mea- 
sure Xmin/(^ ~ 2) so that parameter values coming from 
worse-quality fits caused mainly by the presence of inter- 
lopers were naturally attenuated. 

Fig. [6] illustrates this approach both in the form of the 
final ±3aios(^^) lines on the velocity diagram (left column) 
and the velocity dispersion profiles obtained in subsequent 
steps of this procedure (right column) for halo 6 in three 
projections. The procedure allows to remove on average 61 



8 



R. Wojtak et al.: Interloper treatment in dynamical modelling of galaxy clusters 



0.5 1 1.5 2 
R [Mpc] 



0.5 1 1.5 2 
R [Mpc] 



2 -1 * • V, .'. • * 



0.5 1 1.5 2 
R [Mpc] 



1.2 

I ^ 
^■2 0.8 
o 

=L 0.6 
J 0.4 
0.2 




0.5 1 1.5 2 
R [Mpc] 



1.2 

I ^ 

^f^ 0.8 
b 

n, 0.6 

i 0.4 
0.2 




0.5 1 1.5 2 
R [Mpc] 



1.2 

I ' 

f 0.8 
b 

^ 0.6 

" 0.4 



0.2 



0.5 1 1.5 2 
R [Mpc] 




0=1 
C=5 

c=10 



0.2 0.4 0.6 0.8 1 
R/r^ 

Fig. 7. Boundary lines of the area occupied by halo 
particles with probability — 0.9973 for three values 
of the concentration parameter c = 1,5,10. The velocity 
diagram is assumed to have a cut-off in radius at i?niax = 
Ty and wiim is expressed in units of the local value of the 
projected velocity dispersion trios (i?). 



Fig. 6. Results of the 3a\os{R) method for halo 6 in three 
projections {x, y and z from top to bottom). Left column 
panels show velocity diagrams with final ±3crios(^) lines 
separating interlopers from halo particles. Filled and empty 
circles on velocity diagrams indicate bound and unbound 
particles respectively. Right column panels show the line- 
of-sight velocity dispersion profiles obtained in successive 
steps of the procedure (solid lines from top to bottom). 
The dashed line with error bars is the dispersion profile 
measured for bound particles. 



percent of unbound particles from a given velocity diagram 
with the rate of misidentification of only around 0.2 percent 
(see 3fTios(^) method in Tabled]). Note that fi in this case 
achieves values similar to those obtained with the 3(Tios(10) 
method, while fg is much smaller. 

3.4. Generalized statistical approach 

The methods presented in the above subsection assume im- 
plicitly that the projected distribution of halo particles does 
not depend on the projected radius. In the following, we re- 
formulate this approach properly, taking into account the 
full dependence of the distribution of halo particles on R 
and V. A natural extension of the criterion introduced by 
Yahil & Vidal (1977) is then given by the boundary line 
iviini{R) which determines an area occupied by halo parti- 
cles on the velocity diagram with some probability piim(-R)- 
Conditions for wiim(^) can be written as follows: 



PB.,v[R,v\im{R)]dvdR — CdvdR 



(9) 



where pr.v is the projected probability distribution of halo 
particles and C is its constant value along the boundary line 
itiiim- Second equation in ([9]) is necessary to fully constrain 
the final solution. 



In order to find a useful analytical approximation for 
Pr,v we first assume that the probability of finding a halo 
particle inside an infinitesimal range of radius [i?, R -f di?] 
of a cylinder of observation with radius i?max is given by 



PndR = 1-kR 



Mp{R 



rdR, 



(10) 



where I]m(^, c) and Mp(i?niax,c) are surface density and 
projected mass inferred from this surface density; in our 
case they follow from the NFW density formula (see 
Bartelmann 1996; Lokas & Mamon 2001). To be precise, 
equation pUj) is satisfied best by dark matter particles and 
could be less suitable for description of spatial distribution 
of galaxies. However, cluster data are consistent with galaxy 
distribution being given by the NFW profile and mass-to- 
number density being constant within the virial radius (see 
e.g. Biviano & Girardi 2003) so that formula (|10p seems to 
be useful also in the context of real data. 

Second, we assume that the distribution of the line-of- 
sight velocity at a given radius R (i.e. the conditional prob- 
ability of a particle having a velocity v if it is at projected 
distance R) can be well approximated by a Gaussian dis- 
tribution 



Pvdv = 



1 



27rCTlos(-R, c] 



exp 



2^L(^,c) 



du, 



(11) 



where crios(-R) is the projected isotropic solution to the 
Jeans equation. Although this approximation is not exactly 
valid since departures from Gaussianity are seen both in 
simulated haloes and real clusters (Kazantzidis et al. 2004; 
Wojtak et al. 2005; Hansen et al. 2006; Lokas et al. 2006) 
we claim it is sufficient for our needs. In particular, one can 
show that for isotropic orbits the projected velocity dis- 
tribution following from the distribution function for the 
NFW profile (see Lokas & Mamon 2001; Widrow 2000) is 
remarkably close to a Gaussian in a wide range of radii. 
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Fig. 8. Histograms of the ratio Mp/Myr for different sub- 
samples of particles. Subsamples are identified in the upper 
right corner of each panel: the lower right panel is for all 
particles, upper left for bound particles and the remain- 
ing two for particles inside or 2ry. All histograms are 
normalized to unity. 

Finally, combining the formulae for pn and py we get a 
heuristic expression for the projected probability distribu- 
tion: 



PR.vdRdv = pRPvdRdv. 



(12) 



One can immediately check that the normalization condi- 
tion on the available area of the velocity diagram is satis- 
fied automatically with sufficiently high numerical precision 
since Winax(^) > 4crios(i?): 



R„ 



PRydvdR = 1. 



l.max(i?.) 



(13) 



Note that pn ^ given by (HH), by analogy with results ob- 
tained by Maoz & Bekenstein (1990), maximizes Shannon's 
entropy for known functions Yim{R) and a\os{R)- In the lan- 
guage of the information theory this means that pn.v is the 
most plausible probability distribution given Tim{R) and 

O-los(-R)- 

Substituting into ^ we derive vnmiR) / <Jios{R) 

from the second of these equations. Introducing this ex- 
pression to the first equation we calculate numerically 
the constant C and in the end evaluate the whole profile 
^^iim(^)/o'ios(^)- Fig.[7]shows numerical solutions assuming 
cut-off radius i?max = fy and piun = 0.9973, a value cor- 
responding to the ±3tT range for a Gaussian distribution. 
We see that in general vnm (R) / trios (R) departs significantly 
from the value of 3. Only in the limit of ~ const we re- 
cover the continuous version of the criterion introduced by 
Yahil & Vidal (1977): wiim(i?) = ^<tios{R)- We find that us- 
ing the exact solution to the set of equations ([9]) does not 
improve the performance of the method, the numbers of in- 
terlopers removed in this case are similar as for the simpler 
3a\os schemes (see method in Table [T]). 

3.5. Mass estimators as indicators of interloper fraction 

In this subsection we study the effect of interlopers on the 
values of two standard mass estimators and use the results 




Fig. 9. The ratio of mass estimators Mp and Mvt, their 
derivatives AM / {Arij-cmM) and fractions of unbound fi 
and bound fg particles removed in the procedure based 
on jackknife statistics as a function of the total number of 
removed particles nrom- AH profiles are averaged over the 
whole sample of 30 velocity diagrams and the shaded areas 
indicate the dispersion of values. The solid and dashed lines 
in the lower right panel correspond respectively to Myx and 
Mp. 



to construct another method of interloper removal. The es- 
timators we use are the virial mass Mvt expressed by for- 
mula ([6]) and the projected mass Mp for isotropic orbits 
given by (Heisler et al. 1985; Perea et al. 1990) 



Mf 



32 
nGN 



(14) 



where Vi and Ri are the velocity and projected radius of 
the i-th galaxy, N is the number of galaxies in the sample. 

Both mass estimators are sensitive to the presence of in- 
terlopers, but each of them in a different way (Perea et al. 
1990). Mp is considerably more overestimated than Myr 
so that interlopers effectively give rise to the increase of 
Mp/MvT ratio above 1. In principle, it is impossible to 
relate the value of this ratio to the number of any kind of 
interlopers because of strong degeneracy: a given value of 
Mp/MvT can be reproduced with various numbers of in- 
terlopers and their different distributions on the velocity 
diagram. Nevertheless, it is interesting to study the distri- 
butions of the ratio Mp /Mvt for four different particle sub- 
samples: all particles, particles inside 2ry or ry and bound 
particles. 

The results are shown in Fig.[8]in the form of histograms 
(normalized to unity) constructed from the data in our 30 
velocity diagrams. As we can see, there is no significant 
difference between histograms obtained for all three sub- 
samples with different types of interlopers subtracted. All 
of them have a maximum at around Alp /Mvt = 1 and a 
spread between ~ 0.8 and ~ 1.4. Including unbound parti- 
cles in the analysis, we get a highly asymmetric histogram, 
shifted to the range 1.2 — 2.0 with a maximum at ~ 1.2. 
These results prove that the unbound particles are the ones 
that give rise to the overestimation of the ratio of mass es- 
timators and therefore the ratio Mp/MvT can be a useful 
indicator of the contamination of a sample with unbound 
particles. 
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Fig. 10. Stacked velocity diagrams in projection along x, y and z axes. 



This phenomenological relation between the ratio of 
mass estimators and the presence of unbound particles in 
the sample can be used to construct a procedure which elim- 
inates this kind of interlopers from the velocity diagram. 
The prescription for such interloper identification was pro- 
posed by Perea et al. (1990) and is based on jackknife statis- 
tics. Let {Ri, Vi} he a sequence of data, where i goes from 
1 to n. Following the jackknife technique, we calculate n 
values of both mass estimators which correspond to n sub- 
sequences with one data point excluded. Finally we identify 
as an interloper the particle for which the corresponding 
subsequence is the source of the most discrepant value of 
one of the estimators with respect to the mean value. In 
the next step, the same procedure is applied to a new data 
set with n — 1 particles. 

The main problem of this procedure lies in defining 
properly the condition for stopping the algorithm. In or- 
der to specify it we calculate the ratio Mp/Myr and the 
fractions fi and fg (for interlopers defined as unbound par- 
ticles) as functions of the number of removed particles rircm 
determined by the jackknife technique. Fig. [5] shows the re- 
sults averaged over all 30 velocity diagrams so that the 
shaded areas indicate the dispersion of the values. From 
the behaviour of all profiles we infer that the most ac- 
curate moment of convergence (around nrem/300 = 0.05) 
coincides clearly with a characteristic knee-like point of 
Mp/MvT and IS.M/M profiles. This point is usually more 
recognizable in the case of single velocity diagrams so that 
it can be effectively used to judge where the algorithm 
should be stopped (see Wojtak & Lokas 2006). Taking 
'^rem/300 = 0.05 wc are able to remove on average 65 per- 
cent of unbound particles with around 1 percent of bound 
particles taken for interlopers by mistake (see Mp/Myr 
method in Table [T|). 

4. Indirect methods of interloper treatment 

4.1. The distribution of interlopers 

The key idea of the indirect approach is to treat interlop- 
ers statistically in the proper dynamical analysis so that 
no particles need to be rejected (van der Marel et al. 2000; 
Prada et al. 2003). In all versions of this method one as- 
sumes that the probability distribution p(i?, v) of particles 
seen on the velocity diagram consists of two terms which 
describe the distribution of halo particles and the distribu- 
tion of interlopers. 

Such an analysis requires numerous and regular samples 
of both kinds of particles on the velocity diagram and can- 
not be done for a single object where the distribution of in- 
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Fig. 11. Surface density (left column) and velocity his- 
tograms (right column) of interlopers on stacked velocity 
diagrams in different projections marked in the upper left 
corner of each panel. Solid, dashed and dotted lines corre- 
spond to unbound particles and particles beyond 2ry and r„ 
respectively. Both surface density and velocity histograms 
are normalized to the number of a given type of interlopers 
averaged over 10 velocity diagrams in a given projection. 



terlopers can be highly irregular. We therefore stacked all 30 
velocity diagrams into three composite ones in each projec- 
tion separately. All radii and velocities were rescaled by 
and Vy respectively so that the mass dependence is factored 
out. This procedure is commonly applied to the cluster data 
to improve the statistics and to study the typical proper- 
ties of clusters which are expected to scale with mass (e.g. 
Mahdavi & Geller 2004; Lokas et al. 2006). Fig. [TO] shows 
our three stacked velocity diagrams in projection along the 
x, y and z axes. 

For each of these diagrams, we calculate surface density 
profiles and velocity histograms for all three kinds of inter- 
lopers (Fig. fTTj) . It is clearly seen that unbound particles are 
the type of interlopers with the most uniform surface den- 
sity and velocity distribution, whereas particles from the 
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Fig. 12. The velocity histogram and the fitted probability 
distribution given by p5|) for the outer bin of a stacked 
velocity diagram in projection along the x axis. 



samples beyond 2rv or r^, are considerably concentrated in 
the vicinity of the halo mean velocity (see also Cen 1997, 
Diaferio et al. 1999 and Chen et al. 2006). The uniformity 
of the distribution of interlopers on the velocity diagram 
is the most natural assumption and has been used previ- 
ously in the construction of the probability distribution of 
interlopers (Mahdavi & Geller 2004). We can see that this 
hypothesis agrees well with a distribution of unbound par- 
ticles. 



4.2. Estimation of ttie velocity dispersion profile with a 
uniform background of interlopers 

This method was originally introduced by Prada et al. 
(2003) in the context of measuring the velocity dispersion 
profile of satellite galaxies. Following these authors we as- 
sume that the distribution of particles is given by the sum 
of a Gaussian part for halo particles and a constant back- 
ground describing interlopers. Let fn{v)dv be the condi- 
tional probability of finding any type of particle in the in- 
finitesimal range [v,v + dv] at a given radius R 



fR{v)dv = a{R)p^dv +[1- a{R)] 



dv 

2'^ma, 



(15) 



where is the Gaussian distribution given by (jlip with dis- 
persion equal to aios{R) and the local mean velocity fJ.{R), 
Umax is the maximum velocity available on the velocity di- 
agram and a{R) has a simple interpretation of the proba- 
bility of finding a halo particle at a given radius R. Note 
that in this case the normalization condition is valid only 
for a fixed radius R. All radius-dependent quantities are es- 
timated in radial bins by fitting formula (|15p to a velocity 
histogram in a given bin by minimizing the function. We 
use 10 radial bins with 300 data points in each of them. 

An example of the velocity histogram together with 
the fitted probability distribution given by (1151) is shown 
in Fig. [121 The dispersion profiles obtained in the fitting 
procedure are shown in the left column of Fig. [13] (filled 
circles). For comparison, we also plot velocity dispersion 
profiles of bound particles with shaded areas corresponding 
to ±3(7 departures from the mean. As we can see, all pro- 
files determined in the fitting procedure decline properly 
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Fig. 13. Left column panels show dispersion profiles ob- 
tained by fitting the probability distribution P^]) to ve- 
locity histograms for composite haloes in different radial 
bins (filled circles). Black lines plot dispersion profiles of 
bound particles with shaded strips indicating the 3cr range 
of variability among the velocity diagrams. Right column 
panels show profiles of the probability a{R) measured for 
particles beyond (dotted lines) or 2r^ (dashed lines) and 
for unbound particles (solid lines). Dashed-dotted broken 
lines represent results of the fitting procedure. Each row of 
panels corresponds to a different projection direction of the 
stacked velocity diagram marked in the bottom left corner. 

with radius and trace well the dispersion profiles of bound 
particles. Interestingly, some clearly overestimated values of 
CTios and a probability appear in the same range of radii for 
the case of projection along the z axis which is an effect of 
local irregularities in the velocity distribution of unbound 
particles. 

It seems interesting to compare a{R) profiles obtained 
in the fitting procedure with those directly measured from 
the data for different types of interlopers. Using frequency 
definition of probability a is expressed as 



a{R) = 



Ng{R) 



NgiR) + N,iR)' 



(16) 



where Ng and Ni are numbers of halo particles and in- 
terlopers defined by a given criterion respectively. In the 
right column panels of Fig. [13] we plot a profiles estimated 
in radial bins with interlopers defined as: unbound parti- 
cles (solid lines), particles beyond 2ry (dashed lines) and r„ 
(dotted lines). Profiles obtained in the fitting procedure are 
indicated with broken dashed-dotted lines. We can clearly 
see that, as expected, the fitted a profiles reproduce direct 
measurements best for bound particles. 

We therefore conclude that the favoured group of in- 
terlopers in this approach consists of unbound particles, 
whereas most particles from the outside of the virial re- 
gion, but bound to a given halo contribute to the Gaussian 
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part of the distribution function /_r,(w) and in fact are used 
in the estimation of the dispersion profile. This situation is 
difficult to avoid; wishing to include particles beyond 2r„ 
or ry as interlopers one would have to introduce another 
Gaussian-like distribution of interlopers in velocity space 
(see velocity histograms in Fig. [Tl]). This, however, would 
cause strong degeneracy since both halo particles and in- 
terlopers would be described by very similar distributions. 

4.3. Bayesian technique 

In this subsection we study the elegant approach of sta- 
tistical treatment of interlopers originally proposed by van 
der Marel et al. (2000) and Mahdavi & Geller (2004).' This 
method is based on the Bayes technique which allows us 
to determine the probability distribution in the parame- 
ter space of a particular model given the measured data. 
Consider the parameter set a and data set {a^}. Following 
Bayes theorem the probability of getting certain values of 
parameters a given data sequence {xi} is 



p{a\{xi}) = 



Pja) 



Ilip{xi\a), 



(17) 



where p{a) is the prior on the parameters and p{{xi}) takes 
care of normalization. The combination of p{xi\a) on the 
right-hand side of equation (fTT]) is the likelihood, while 
p(a|{a^}) is the posterior probability. Obviously, each of the 
probability distributions introduced above is normalized to 
unity in the available part of the corresponding space 



-oo p+oo 



/-too 
p{a\{xi})dai...dan = 1 
-oo 

/+00 r+oo 
... / p{xi\a)dxi...dxn = 1- 
-oo J ~oo 



(18) 
(19) 



Van der Marel et al. (2000) proposed to restrict the con- 
siderations to the velocity space {xi = Vi). The probabil- 
ity p{vi\a) was given by (|15p with a more detailed formula 
for py and the assumption that a{R) « const, which was 
found to be roughly consistent with the data. However, as 
we have seen in Fig. [13] the dependence of the probability 
a on radius is in general not negligible. The most natural 
way to take this fact into account is simply to consider the 
probability on the whole projected space with Xi — {Ri, Vi) 



p{Ri,Vi\{ad,ap))dRdv = apf{Ri,Vi\ad)dRdv 

/ ^ R 

+ (l-ap) 



(20) 



^max'^^max 



dRdv, 



where /(-R, v) is the projected distribution function of halo 
particles, is a set of dynamical parameters and ap is an 
additional free parameter describing the probability that 
a particle found at any radius and with any velocity is a 
halo particle. The last term in (PD|) describes the uniform 
distribution of interlopers both in position on the sky and 
in velocity space. 

We performed the analysis of our three stacked haloes 
using the probability distribution ([20]) with f{Ri,Vi\a) 
given for simplicity by the formula (|12[) parametrized by 
the c parameter. Since we are dealing with stacked haloes 
and we assume isotropic orbits (as shown to be the case for 
early-type galaxies in Coma and other clusters, see Lokas 
& Mamon 2003; Biviano & Katgert 2004), the only free 
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Fig. 14. Left column panels show the results of the 
Bayesian analysis in the form of contours in the c-a^ pa- 
rameter plane corresponding to 68.3, 95.4 and 99.73 percent 
confidence levels. In the right column we plot a{R) prob- 
ability profiles calculated for different types of interlopers: 
particles beyond (dotted lines) or 2r« (dashed lines) and 
unbound particles (solid lines). Shaded areas indicate a{R) 
values corresponding to c and ap parameters from the in- 
side of the 99.73 percent probability contours. The rows of 
panels from top to bottom are for the x, y and z projections 
respectively. 



parameters in this analysis are the concentration and the 
probability ap. Results of this analysis are shown in the left 
column of Fig. [TJ] in the form of contours corresponding to 
68.3, 95.4 and 99.73 percent confidence levels drawn in the 
c-Op parameter plane. 

As we can see from Fig. [141 the preferred values of c 
lie slightly below the value of the concentration parameter 
(estimated from the full 3D data) for the composite cluster, 
c = 6.9, which is probably due to our simplistic probability 
distribution f{R,v). However, we believe that the results 
concerning Up illustrate well the general features of inter- 
loper treatment in this approach since the distribution of 
interlopers used in formula (1^0]) is independent of the dy- 
namical model. 

It is interesting to compare Up with abundances of dif- 
ferent types of particles seen on velocity diagrams and the 
mean value of a obtained in the previous subsection. All 
these values are listed in Table [21 where F^y^^^ and -F<2r„ 
are the abundances of bound particles and particles within 
2r„ respectively. Comparing these results, we notice an ex- 
cellent agreement between ap and a. This fact is under- 
standable since a was estimated in radial bins with the 
same number of particles inside each of them. Consequently 
one expects that the mean value of a measures the proba- 
bility that a particle randomly chosen from the velocity dia- 
gram is a halo particle which is exactly the same probability 
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Table 2. Probabilities (ctp and a) of finding a halo parti- 
cle on the velocity diagram estimated via indirect methods 
of interloper treatment. F^y^^^ and -F<2r„ are fractions of 
bound particles and particles within 2ry respectively. 
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as described by ap. We also find that ap is approximately 
equal to the abundance of bound particles. This means that 
in this approach mainly unbound particles contribute to the 
distribution of interlopers. Some departure from this rule 
can be seen for the z-axis projection. However, if we take 
into account that the velocity distribution of unbound par- 
ticles in this case is slightly peaked close to u ~ (see 
the bottom right panel of Fig.[Tl|) the situation is clarified: 
some of these interlopers are identified as halo particles in- 
creasing ap. 

Although we used here only the global probability ap of 
finding a halo particle in the velocity diagram, it is easy to 
reproduce the a{R) profile. In this case equation (fTB|) has 
its counterpart in the form 



a{R) 



(21) 



where = const is the surface density of interlopers re 
lated to Up by 



27ri?EM(i?, c)dR 



" 2TTRT.MiR,c)dR 



(22) 



Therefore given ap and c one is able to calculate a{R) pro- 
file. The results are plotted in the right column panels of 
Fig. [T3] in the form of shaded areas which indicate a{R) 
values calculated for c and ap parameters from the inside 
of the probability contour corresponding to 99.73 percent 
confidence level. We can clearly see that the area available 
for a{R) includes a profile calculated for interlopers as un- 
bound particles. 

We therefore confirm the results of the previous subsec- 
tion that the uniform distribution of interlopers in the ve- 
locity diagram is mainly reproduced by unbound particles, 
whereas most of bound particles from the outside of the 
virial sphere contribute to the distribution function of halo 
particles. We emphasize that only indirect methods are able 
to take into account the presence of low velocity interlopers 
in the central part of the velocity diagram. Consequently, 
they can potentially include all unbound particles in the 
interloper background, which is beyond the reach of any 
direct scheme. For instance, applying 'ymax(l) method, the 
most successful among direct ones, to the composite ve- 
locity diagrams we found that 21, 16 and 44 percent of 
unbound particles (for x, y and z projection respectively) 
are not identified and remain in the final samples of halo 
members resulting in the contamination on the level of 0.9, 
1.7 and 4.4 percent respectively. 



5. Discussion 

We have studied different approaches to the treatment of 
interlopers in the analysis of kinematic data for galaxy clus- 
ters. For the direct methods of interloper removal their ef- 
ficiency was measured by simple parameters: the fraction 
of removed interlopers /i, the fraction of removed mem- 
bers fg and the fraction of non-identified interlopers re- 
maining in the final samples fh- The values of these pa- 
rameters obtained by avergaing over 30 velocity diagrams 
studied here are given in Table [TJ We can see that the 
highest fi (73 percent for unbound particles) is reached 
by i'max(l) method originally proposed by den Hartog & 
Katgert (1996) although all the other methods considered 
(except for Wmax(2)) have fi on the level of 60 percent. The 
differences in mean fi between the methods are small com- 
pared to the dispersion obtained by averaging the results 
over 30 velocity diagrams which is about 20-30 percent in 
the case of unbound particles. In all direct methods of in- 
terloper removal fi equals approximately 40 percent and 20 
percent respectively for particles beyond 2r^ and r„ so the 
methods are not efficient in removing them. This occurs 
because particles from the close surroundings of the virial 
region are significantly concentrated around the halo mean 
velocity and may in fact be good tracers of the halo po- 
tential. All the methods show rather low fractions of mem- 
bers misidentified as interlopers (below 2 percent for all 
types of interlopers with a dispersion of about 2 percent). 
Furthermore, final samples of halo members include on av- 
erage 2-4 percent of unbound particles with dispersion of 
around 5 percent. 

One may wonder if this is the most reliable way to com- 
pare the different methods. Since our purpose is to maxi- 
mize fi and minimize fg and fh one could construct some 
combination of these parameters assigning them different 
weights; we have not done this in order not to obscure the 
picture. One can also ask whether removing 60 percent of 
interlopers by one method cannot be better than remov- 
ing 70 percent by another when it comes to the final esti- 
mation of the cluster parameters. It could happen that a 
smaller number of removed interlopers could lead to bet- 
ter estimates of mass because the interlopers were actually 
those causing the largest bias while larger number of inter- 
lopers can in principle lead to more biased results because 
the less significant interlopers were removed. Fortunately, 
this is not the case, the interlopers removed by the differ- 
ent methods are mostly the same, the difference is usually 
in the border-line particles which do not significantly con- 
tribute to the bias. We have also verified that the final, 
cleaned samples depend very weakly on the initial cut-off 
in velocity: reducing the cut-off to 3000 km s~^ from 4000 
km with respect to the cluster mean produces almost 
identical final samples differing by only 1-2 particles in a 
few out of 30 velocity diagrams. 

Apart from the quantitative measures presented in 
Table [T] the choice between different methods of interloper 
removal is largely a matter of subjective preference. One 
could try to judge the methods by how strong their un- 
derlying assumptions are but the methods rely on such 
different assumptions that it is difficult to compare them. 
Another possibility is to look at their convergence. Here 
again Wmax(l) method is recommended: it always converges, 
contrary e.g. to second most effective method Mp/Myr 
based on the ratio of mass estimators where the procedure 
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has to be stopped at a rather arbitrary point. Furthermore, 
as shown by van Haarlem et al. (1997), Vmax(l) approach 
leads to the data sample which reproduces real velocity dis- 
persion better than the one obtained after applying 3(Tios 
algorithm in the original form proposed by Yahil & Vidal 
(1977). It is also worth noting that t'max(l) method does 
not involve any parameters or characteristic scales which 
could restrict its application to cluster-size objects. It has 
been recently successfully applied to tidally stripped dwarf 
spheroidal galaxies where it allows to clean the stellar kine- 
matic samples from interlopers originating from the tidal 
tails and the Milky Way (Klimcntowski et al. 2006). 

An ultimate verification of the methods would of course 
come from the dynamical modelling performed on the 
cleaned samples. Here, however, other issues come into play: 
the quality of the final results of the modelling depends 
also on the actual method used and particular properties 
of an object (whether it is in equilibrium, whether it is 
spherically symmetric, how well is the kinematics sampled). 
These sources of error may have more impact on the final 
result than the differences between the methods of inter- 
loper removal studied here (see the recent study by Biviano 
et al. (2006) who have analyzed, after interloper removal 
with ?;max(l) method, different observational effects affect- 
ing the determination of cluster mass and velocity disper- 
sion). However, if interlopers were not removed at all the 
result would be very strongly biased and the contamination 
could become the main source of error. In any case it would 
be very difficult to disentangle the effect of interlopers from 
other sources of uncertainties. 

Comparison between many modelling approaches pos- 
sible is beyond the scope of this paper and it was not our 
purpose to provide here; such a final answer. For the study of 
the performance of one of the dynamical analysis methods 
based on fitting the velocity dispersion and kurtosis pro- 
files for clusters we refer the reader to Sanchis et al. (2004) 
and Lokas et al. (2006). Examples of such full dynamical 
modelling, including the interloper removal as an impor- 
tant first step, in the case of Abell 576 and other clusters 
with significant background contamination can be found in 
Wojtak & Lokas (2006). 

For the purpose of the present study we performed a 
simple test taking the cleaned particle samples for our 30 
velocity diagrams obtained with Wmax(l) method and fit- 
ting the velocity dispersion profiles to the solutions of the 
Jeans equation assuming isotropic orbits and estimating 
the virial masses and concentrations. Averaging over 30 
diagrams we find that the ratio of the estimated virial 
mass to the real one measured from the 3D information 
is My/ My ^true = 0.86 ± 0.23 while for concentrations we 
get c/ctrue = 1-76 ± 1.01. Interestingly, for the samples of 
bound particles the results are very similar: My/My^tiuc = 
0.85 ± 0.18, c/ctruo = 1-91 ± 1.31. These values should be 
compared to the parameters obtained for the samples of 
all particles (without application of any interloper removal 
scheme): M„/M„,true = 1-64 ± 0.92, c/ctrue = 0.34 ± 0.31. 
We therefore conclude that using the samples cleaned with 
fmax(l) method is equivalent to working with only bound 
particles. The significant bias still present in both cases, 
especially for concentration, can be explained by depar- 
ture from isotropy of particle orbits in the simulated haloes 
which should be taken into account in more complete mod- 
elling. 



We have also considered indirect methods of interloper 
treatment where their presence is accounted for in a sta- 
tistical way. Here we have found that, contrary to the as- 
sumptions of van der Marel et al. (2001), the probability of 
a given galaxy being an interloper is not independent of the 
projected radius but increases with it. We have also veri- 
fied the applicability of the approach of fitting a Gaussian 
plus a constant to the velocity distribution in galaxy clus- 
ters as a method to account for interlopers in estimating 
the velocity dispersion profile originally proposed by Prada 
et al. (2003) to study the velocity distribution of satellites 
around giant galaxies. 

The main disadvantage of the indirect methods is that in 
order to reliably estimate the parameters they require large 
kinematic samples which can only be obtained by stacking 
data coming from many objects. This procedure, although 
commonly used, can be dangerous, because it is not clear 
how the distances and velocities should be scaled. While 
in the case of distances the choice of the virial radius as 
the scaling parameter seems rather obvious, in the case of 
velocities it is less clear whether one should use velocity 
dispersion, circular velocity at the virial radius or perhaps 
the maximum circular velocity and still all would be sub- 
ject to imcertainties due to modelling of single clusters. In 
addition, the underlying assumption of indirect methods is 
that the background of interlopers is uniform which is never 
exactly the case due to clustering. 
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